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1 Introduction 



• The lattice Boltzmann scheme is a numerical method for simulation of a wide family of 
partial differential equations associated with conservation laws of physics. The principle is 
to mimic at a discrete level the dynamics of the Boltzmann equation. In this paradigm, the 
number f(x, t) dx dv of particles at position x, time t and velocity v with an uncertainty 
of dx dv follows the Boltzmann partial differential equation in the phase space (see e.g. 
Chapman and Cowling [7]): 

(1) % + v.VJ = Q(f) . 



• Note that the left hand side is a simple advection equation whose solution is trivial 
through the method of characteristics: 

(2) f(x,v,t) = f{x-vt,v,0) if Q(/) = 0. 

Remark also that the right hand side is a collision operator, local in space and integral 
relative to velocities: 

(3) Q{f){x, v,t) = J C(f(x, w, t),x, v, t) dw, 

where C(») describes collisions at a microscopic level. Due to microscopic conservation of 
mass, momentum and energy, an equilibrium distribution f eq (x, v, t) satisfy the nullity 
of first moments of the distribution of collisions: 



Q(D(x,v,t) I v | dv = 0. 

Such an equilibrium distribution / eq satisfies classically the Maxwell-Boltzmann distri- 
bution. 

• The lattice Boltzmann method follows all these physical recommandations with spe- 
cific additional options. First, space x is supposed to live in a lattice C included in 
Euclidian space of dimension d. Second, velocity belongs to a finite set V composed by 
given velocities Vj (0 < j < J) chosen in such a way that 

x G C and Vj G V =^ x + At v j G C , 

where At is the time step of the numerical method. Then the distribution of particles, 
/, is denoted by fj(x, t) with < j < J, x in the lattice £ and t an integer multiple 
of time step At. 
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• In the pioneering work of cellular automata introduced by Hardy, Pomeau and De 
Pazzis [24], Frisch, Hasslacher and Pomeau [18] and developed by d'Humieres, Lallemand 
and Frisch [IB], the distribution fj(x, t) was chosen as Boolean. Since the so-called 
lattice Boltzmann equation of Mac Namara and Zanetti [35], Higuera, Succi and Benzi 
[27] . Chen, Chen and Matthaeus (HJ, Higuera and Jimenez [26] (see also Chen and Doolen 
[9]), the distribution •) takes real values in a continuum and the collision process 
follows a linearized approach of Bhatnagar, Gross and Krook [4]. With Qian, d'Humieres 
and Lallemand [38], the equilibrium distribution / eq is determined with a polynomial in 
velocity. In the work of Karlin et al [29], the equilibrium state is obtained with a general 
methodology of entropy minimization. 

• The numerical scheme is defined through the evolution of a population fj(x, t), with 
x E £ and < j < J towards a distribution fj(x, t + At) at a new discrete time. 
The scheme is composed by two steps that take into account successively the left and 
right hand sides of the Boltzmann equation (pQ). The first step describes the relaxation 
/ — > f* of particle distribution / towards the equilibrium. It is local in space and 
nonlinear in general. D. d'Humieres first introduced in [ 1 1 J the fundamental notion of 
moments in the context of lattice Boltzmann schemes. He defines an invertible matrix 
M with (J+ 1) lines and (J + l) columns and the moments m through a simple linear 
relation 

j 

(4) m h = Y,M kj fj, 0<k<J. 

3=0 

• The first N moments are supposed to be at equilibrium: 

(5) m* = m i = mf l = Wi, 0<i<N-l 

and we introduce the vector W G M. N of conserved variables composed of the Wi for 
< % < N — 1: Wi = m^, < % < iV — 1. The first moments at equilibrium are 
respectively the total density 

j 

(6) P = E^> 

3=0 

momentum 

j 

(7) q a = J2 v ?f3, l<a<d 

3=0 

and possibly the energy for Navier-Stokes fluid simulations. In consequence, we have 

(8) M 0j = 1, 0<j<J 

(9) M aj = v J , 1 < a < d , < j < J . 
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For the other moments, we suppose given ( J + 1 — N) (nonlinear) functions 

(10) R N 3Wi — > G k (W) eR, N<k<J 
that define equilibrium moments m e k q according to the relation 

(11) m e k q = G k (W), N < k < J . 

Note also that more complicated models have been developed in Yeomans's group (see 
e.g. Marenduzzo at al [36]) for modelling of liquid crystals. 

• The relaxation process is related to the linearized collision operator introduced at 
relation (J3]). In particular intermolecular interactions (Maxwell molecules with a 1/r 4 
potential), the collision operator is exactly solvable in terms of so-called Sonine polyno- 
mials (see e.g. Chapman and Cowling |7J) and the eigenvectors are known. Moreover, 
the discrete model is highly constrained by symmetry and exchanges of coordinates. In 
the work of d'Humieres [TT], relaxation parameters (also named as s-parameters in the 
following) Sk (N < k < J) are introduced, satisfying for stability constraints (see e.g. 
|30j) the conditions 

< s fc < 2 , N <k< J . 
Then the nonconserved momenta m* k after collision are supposed to satisfy 

(12) m* k = m k + s k (m e k q - m k ) , k> N 

and we will denote by S the diagonal matrix of order J+l—N whose diagonal coefficients 
are equal to s k : 

(13) Ski = 5 ke si>, k, £> N 

with Ski the Kroneker symbol equal to 1 if k — I and null in the other cases. Remark that 
this framework is general: when the matrix S is proportional to identity, the d'Humieres 
scheme degenerates to the popular "BGK" method characterized by a "Single Relaxation 
Rate". In this particular case the relaxation operator is diagonal and there is no particular 
diagonalization basis to work with. The distribution /* after collision is reconstructed 
by inversion of relation (0]): 

J 

(14) /; = ^M-/m £ *, 0<j<J. 

• We suppose also that the set of velocities V is invariant by space reflection: 

Vj e V 3 £ e {0, . . . , J}, Vi = -Vj , v e e V . 
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The second step is the advection that mimic at the discrete level the free evolution through 
characteristics ((2j): 

(15) f j {x,t + At) = f*{x-v j At,t), x e C, < j < J, Vj e V. 

Note that all physical relaxation processes are described in space of moments. Neverthe- 
less, evolution equation ffl~5"l) is the key issue of forthcoming expansions. 

• The asymptotic analysis of cellular automata (see e.g. Henon [25]) provides evidence 
supporting asymptotic partial differential equations and viscosity coefficients related to 
the induced parameter defined by 

11 

(16) o k = - - - . 

The lattice Boltzmann scheme (111) to (TT5l) has been analyzed by d'Humieres [UJ with a 
Chapman-Enskog method coming from statistical physics. Remark that the extension of 
the discrete Chapman-Enskog expansion to higher order already exists (Qian-Zhou [39] . 
d'Humieres [12]). But the calculation in the nonthermal case (N > 1) is quite delicate from 
an algebraic point of view and introduces noncommutative formal operators. Recently, 
Junk and Rheinlander [28] developed a Hilbert type expansion for the analysis of lattice 
Boltzmann schemes at high order of accuracy. We have proposed in previous works [T4l[T5] 
the Taylor expansion method which is an extension to the lattice Boltzmann scheme of 
the so-called equivalent partial differential equation method proposed independently by 
Lerat and Peyret [33] and by Warming and Hyett [48J. In this framework, the parameter 
At is considered as the only infinitesimal variable and we introduce a constant velocity 
ratio A between space step and time step: 

v ' At 

The lattice Boltzmann scheme is classically considered as second-order accurate (see e.g. 
[30]). In fact, the viscosity coefficients \i relative to second-order terms are recovered 
according to a relation of the type 

fx = (\ 2 Atak 

for a particular value of label k. The coefficient ( is equal to | for the simplest models 
that are considered hereafter. 



• A natural question is to extend this accuracy to third or higher orders. In the case 
of single relaxation times (the BGK variant of d'Humieres scheme), progresses in this 
direction have been proposed by Shan et al [1H 25] and Philippi et al [37] using Hermite 
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polynomial methodology for the approximation of the Boltzmann equation. The price to 
pay is an extension of the stencil of the numerical scheme and the practical associated 
problems for the numerical treatment of boundary conditions. Note also the work of 
the Italian team (Sbragaglia et al [42], Falcucci et al [17]) on application to multiphase 
flows. In the context of scheme with multiple relaxation times, Ginzburg, Verhaeghe 
and d'Humieres have analyzed with the Chapman-Enskog method the "Two Relaxation 
Times" version of the scheme [22], [23]. A nonlinear extension of this scheme, the so-called 
"cascaded lattice Boltzmann method" has been proposed by Geier et al [19]. It gives also 
high order accuracy and the analysis is under development (see e.g. Asinari [3]). The 
general nonlinear extension of the Taylor expansion method to third-order of accuracy 
of d'Humieres scheme is presented in [16]. It provides evidence of the importance of the 
so-called tensor of momentum- velocity defined by 

j 

(18) A{ p = J2 M kjM PJ Mr/, 0<k,p,£<J. 

3=0 

Moreover, it shows also that for athermal Navier Stokes equations, the mass conservation 
equation contains a remaining term of third-order accuracy that cannot be set to zero by 
fitting relaxation parameters |16| . 

• Our motivation in this contribution is to show that it is possible to extend the order of 
accuracy of an existing a priori second-order accurate lattice Boltzmann scheme to higher 
orders. We use the Taylor expansion method [15] to determine the equivalent partial 
differential equation of the numerical scheme to higher orders of accuracy. Nevertheless, 
it is quite impossible to determine explicity the entire expansion in all generality in the 
nonlinear case. In consequence, we restrict here to a first step. We propose in the following 
a general methodology for deriving the equivalent equation of the d'Humieres scheme at 
an arbitrary order when the collision process defined by the functions Gk of relation 
(TTOl) are linear. This calculation leads to explicit developments that can be expanded 
with the help of formal calculation. This work is detailed in Section 2. In Section 3 
we apply the general methodology to classical linear models of thermics and linearized 
athermal Navier Stokes equations. We treat fundamental examples from one to three space 
dimensions. When it is possible, the equivalent partial equivalent equations are explicited. 
In Section 4, we use the fourth-order equivalent equation of two and three-dimensional 
models to enforce accuracy through a proper choice of "quartic" parameters. For a scalar 
heat equation, the effect of the precision of the numerical computation of eigenmodes is 
presented. For linearized athermal Navier Stokes equations, we propose a method for 
enforcing the precision of the eigenmodes of the associated partial differential equation. 
First numerical results show that for appropriate tuning values of the parameters, fourth- 
order precision is achieved. 
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2 Formal development of linearized d'Humieres scheme 

• In what follows, we suppose that the collision process is linear i.e. that the G k 
functions introduced in (fTOl) (fTTi) are linearized around some reference state. With this 
hypothesis, we can write: 

N N 

G k {W) = ^G^Wj = J2 G V m i> k ^ N - 

To be precise, putting together relations (fTTl) and (fT2l) . there exists a (J + 1) x (J + 1) 
matrix \1/ such that the collisioned momentum m* defined in (TT2]) is a linear combination 
of the moments before collision: 

J 

(19) m* = \l/»m, m* k = ^J^W m i ■ 

Of course, the conservation ([9]) implies that \I/ has a structure of the type 



(20) * 



I 



The top left block of the right hand side of (1201) is the identity matrix of dimension N 
and the bottom left block is described through the G k functions introduced in (fTOl (TTTT) : 

(21) ee m kj = s k G kj , j < N , k > N . 



The bottom right block of the right hand side of (1201) contains the coefficients 1 — s k 
(k > N) related to relaxation (TT3l) . 

• In order to make our result explicit, we need some notations. We introduce multi- 
indices 7, 5, e in {1, . . . d} q in order to represent multiple differentiation with respect to 
space. If 

«i times ay times 

then 

Oy EE 



dx* 1 dxf 
and we denote by |7| the length of multi-index 7: 

1 7 j EE «! + ••• + a d . 

Then thanks to the binomial formula for iterated differentiation, we can introduce coeffi- 
cients in order to satisfy the identity 

d 

(22) (-J2 M ^ d °) q = H P ^ 9 t 

a=l \y\=q 
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for any integer q. 

• We first establish that at first-order of accuracy, we have a representation of noncon- 
served moments in terms of conservative variables: 

j 

(23) m k = J2 B °kj W J + °( At ) > k > N . 

3=0 

with 

(24) B° kj = -^ kj , k>N,0<j<N-l. 

Sk 



We have also the first-order conservation law 



(25) ^ + XX W = O(At), 0<z<iV-l 



l7l=l 

with coefficients AJj given according to 

J -, 

(26) A% s^^^ + E^r^)' l7l=l, 0<i,j<JV-l. 

p=0 £>Af ^ 

The proof of this result and those that follow of this Section are detailed in Appendix A. 

• The expansion of moments (1231) can be extended to second-order accuracy: 

(27) m k = Atbl B l d i W J + °( At2 ) • 

o<H<i 

with 

( 1 1 3 1 J 

(28 , 

[ |7|=1, fc>JV, 0<j<iV-l. 

• Then we extend the previous expansions (1251) and fl27|) to any order a. By induction, 
we establish that we have an equivalent partial differential equation of the form 

N-l 

- 4- V AfH- 1 

dt 



(29) ^ + Yl At ' 7hl E ^ = °( Ata ) > < z < AT - 1 

1<|t|<ct j=0 



and an expansion of nonconserved moments as 

2V-1 

(30) m k = At ' 7 ' E B ki d f W i + °( A ^ +1 ) > k> N . 

0<| 7 |<a j=0 
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with the following recurrence relations for defining the coefficients Aj- and Bjj- 



(31) 
(32) 



Clf = AJ j7 0<i,j<N-l, 

- E fc?f^, 

6>q,e>l,5+E=i £=0 

2 < q+ 1 < M , < Z, j <N-l 



(33) 



(34) 



(35) 



l7l 1 



g=2 



E 



J J J J 

E E E nrn Mrf M * *v p " ^ ' 



l<|5|<l7l>0<|e|<|7|-l,(5+e=7 ^=0 p=0 r=0 1 1 

& > A , < j < N - 1 , 
j 



r)0,7 _ n7 



n9 +1 '7 



E E «tf ^ 

|5|>9, |e|>1,<5+£=7 ^=0 

i<9 + i<ItI, fc>iv, 0<j<JV-l, 



(36) ^ 



^ ~ Sfe ( E q \ D W 

J J J v 

EEE t^^m^p^O 



1<9<|7| 

+ E 

1<I*I<|7|,0<|e|<|7|-1,5+£=7 fcO p=0 r=0 



& > A , < j < A - 1 



• Note that the results (1331) and (1361) are coupled through the relations (13T1) (1321) (1341) 
and (1351) . For example, the evaluation of coefficient ' 7 uses explicitly , the eval- 
uation of A?- uses and the computation of B^ is impossible if .D^ 7 is not known. 
The proof is detailed in Appendix A. It is an elementary and relatively lengthy algebraic 
calculation. In particular, our mathematical framework is classical: all differential opera- 
tors commute and the technical difficulties of noncommutative time derivative operators 
associated with the use of formal Chapman-Enskog method [12] vanish. As a result, the 
general expansion of a linearized d'Humieres scheme at an arbitrary order can be obtained 
by making explicit the coefficients AJ- and B\-. Remark that the hypothesis of linearity 
allows making the above formulae explicit and we have done this work with the help of 
formal calculation. Nevertheless, it is always possible to suppose that the Gk functions 
are linearized expansions of a nonlinear equilibrium. In this case, the previous equivalent 
high order partial differential equations (1291) give a very good information concerning the 
behavior of the scheme. 
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3 Equivalent Thermics and Fluid equations 

• We make explicit in this section the fourth-order equivalent equation fl29l) of some lat- 
tice Boltzmann schemes for two fundamental problems of mathematical physics: thermics 
and linearized athermal Navier-Stokes equations. We treat first advective thermics in one 
space dimension with the so-called D1Q3 lattice Boltzmann scheme. In order to obtain 
results presentable on a sheet of paper, we simplify the model and omit the advective 
term for two (D2Q5) and three (D3Q7) space dimensions. Secondly we study linearized 
athermal Navier-Stokes equations in one (D1Q3), two (D2Q9) and three (D3Q19) space 
dimensions. Note that we have to define precisely our results. First the numbering of 
degrees of freedom via corresponding graphics is specified; see Appendix B. The choice 
of moments, id est the M matrix, is also made precise in Appendix B. Secondly the $ 
matrix of relation (fT9l) is specified, later in this section. 




• D1Q3 for advective thermics at fourth-order 

For a thermics problem, we have only one conserved quantity. Then A = 1 in relation ([5|). 
The two nonconserved moments (momentum q eq and energy e eq ; see ([781) ) at equilibrium 
are supposed to be linear functions of the conserved momentum p: 

(37) q^ = u\p, e ^ = a^p. 

Due to fl2T]) and (137]) . the matrix ^ for dynamics relation (fT9l) is given according to 

/ 1 

* = siuX 1 - Si 
\a s 2 A 2 /2 

We determine without difficulty the equivalent partial differential equation for this lattice 
Boltzmann scheme at order four, to fix the ideas. For i — 1, 2, we introduce cr, from 
relaxation time Sj according to relation (fT6|) . When a drift in velocity u is present, note 
that the diffusion coefficient is a function of mean value velocity. We have 

~dt ~dx ~ ai ( a - U ^ + K *—W +K *—d^ = °( A ^ 

with parameters and k 4 given according to 

k 3 = -u (2 {I- 12 al)u 2 + l-3a - I2a x a 2 {l-a) + 24 a? a) 

k 4 = (-9 + 60a^)(Tin 4 + (-5(l-3a)(j 1 - 3(l-a)a 2 + 
+ 12 (1 - a) 01 a\ + 36 (1 - a) o\ a 2 - 72 af a) u 2 
+ a ai (2 - 3 a - 12 (1 - a) a x a 2 + 12 a af) . 
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If u = 0, then k 3 = and the scheme is equivalent to an advection-diffusion equation 
up to third-order accuracy. In this particular case, the scheme is fourth-order accurate in 
the previous sense if we set 

2 -3a + 12 a a\ 



°2 



12 ^{1 -a) 



• D2Q5 for pure thermics at fourth-order 

We have J = 4 and N = 1. The equilibrium energy (momentum m 3 in (1791) with the 
labelling conventions of Section 1) is the only one to be non equal to zero. The matrix 1 !/ 
of relation (fT9l is now given by the relation 



(39) 



/ 1 











\ 






























as 3 








l-s 3 





V 











1-sJ 



We have developed the conservation law up to fourth-order: 

dp A 2 At 
dt ' 



01 (4 + a) 



(40) 
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At 3 A 4 , A 
+ cti (4 + a 



d 2 p d 2 p 



dx 2 

K 40 



d 4 p d 4 p 



+ K 22 



dx 2 dy r " 



1200 ' A v " ' "' ^dx A 
and the k coefficients are explicited as follows: 

(41) n A0 = 8 - 3a + 12 (a + 4) of - 12(1 -a)a x a z - 60aia 4 

(42) k 22 = -6 (a + 4) + 24 (a + 4) a 2 - 24 (1 - a) a x a 3 + 120 a t a 4 



0(At 4 



• D2Q9 for advective thermics at fourth-order 

The lattice Boltzmann model D2Q9 for a passive scalar (see Chen, Ohashi and Akiyama 
[TO] . Shan [43], Ginzburg [20]) is obtained from the D2Q5 model by adding four velocities 
along the diagonals (Figurel2Tl right). The evaluation of matrix M is absolutely nontrivial 
and is precised at (IBTil) . The dynamics is given by 



(43) 



/ 1 























\ 


u A si 


Isi 























v\si 

























A3 «3 








1-S3 

















(24 S4 











1— S4 














a 5 us 5 














1-S5 











a & vs b 

















1-S5 








a 7 Sj 




















1-S7 





\ a 8 s 8 























ls 8 ) 
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The coefficients 03 to as in relation (1431) are chosen in order to obtain the advection 
diffusion equation at order 2: 

We have precisely: 

a 3 = 3 (w 2 + v 2 ) — 4 + 6 £ , 07 = u 2 — v 2 , a 8 = uv 
as explained in [16]. When u — v — 0, the equation (JUT) takes the form 

I " A2 « At (0 + 0) + T ( K « (0 + 0) + A) = ° (At4) 

with coefficients k 4 o and k 2 2 evaluated according to 

k 40 = (J\ ha 5 (0-7 - 0-3) (a 4 - 4) + 6£ (l - a x a 7 - 5axa 3 + 2a 5 (a 7 - cr 3 )) ) 

K 22 = 2 (ai + a 5 - 2 cri cr 5 (cr 3 + cr 7 + 4 cr 8 )) (a 4 - 4) 

+12 £ (o- 5 + 3 cxi - 2 ai a 5 (a 3 + a 7 ) - 2 a t a 3 a 5 -8 o x cr$ {a x + cr 5 ) + o\ a 7 ) . 

Remark that the equivalent partial differential equation of this general lattice Boltzmann 

scheme has been exactly derived in a complex case where all the time relaxations are a 

priori distinct. The coefficients K40 and K22 of the fourth-order terms are polynomials 

of degree 3 in the a coefficients. When we make the "BGK hypothesis" id est that all the 

cr coefficients are equal, a first possibility for killing the coefficients k 40 and k 22 is given 

by: ^ 

o~\ = Ox = (T 3 = (T 4 = <7 5 = 07 = er 8 = - £ = 0. 

6 

We observe that this choice of parameters is without any practical interest because the 
diffusion term in fT44l) is null. We observe that a second possibility 

_ 2 1-6(7? _ n l-2a 2 

4 " 3 l-8a 2 ' ° 4 " 2 1-8(7? 

induces also a fourth-order accurate lattice Boltzmann scheme. If we replace the strong 
"BGK hypothesis" by the weaker one associated to "Two Relaxation Times" as suggested 
by Ginzburg, Verhaeghe and d'Humieres in [22] [23], id est 

0~\ — 05) 03 = cr 4 = a 7 = cr 8 , 

we can achieve formal fourth-order accuracy for 
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• D3Q7 for pure thermics 

For three-dimensional thermics, one only needs a seven point scheme and use the so-called 
D3Q7 lattice Boltzmann scheme whose stencil is described in Figure [22l The matrix M is 
given at relation (18T1) . The dynamics of this Boltzmann scheme uses the following matrix 
for computation of out of equilibrium moments, according to relation (Tl~9l : 

\ 




l-si 






/ 1 














l-s 1 











l-si 












1 — s 4 











1 — S 4 














56/ 



The equivalent thermal scalar conservation law now takes the following form at fourth- 
order of accuracy: 

Qp \2 A+ A+3 \4 / ,<}4, 



A 2 At . . At 3 A 4 , 

dt -^r^^ + ^ A P + Tm ai{a + 6) 

<9 4 p 

+ K 220 ' 

where the k coefficients are given by 

(45) k 400 = 8 - a + 4 o\ (a + 6) 

(46) k 220 = - 2 (a + 6) + 8 o\ (a 



(dp 



d p d 



d 4 p 

dx 2 dy 2 ' dy 2 dz 



+ 



f _&P_\ 

dz 2 dx 2 ) 



0(At 4 



- 56 o"i cr 4 — 4 (1 — a) U\ a% 
6) + 56 <J\ a 4 — 8(1 — a) o x ■ 



• After these examples where only one partial differential equation is present, we con- 
sider the case of two (D1Q3), three (D2Q9) or four (D3Q19) partial differential equations 
"emerging" from the lattice Boltzmann algorithm. These equations model macroscopic 
conservation of mass and momentum of a linearized fluid in our approach in this contri- 
bution. 

• D1Q3 for athermal linearized Navier-Stokes at fifth-order 

We have in this case two conservation laws (N = 2 in relation and the equilibrium 
energy is supposed to be given simply by 



(47) 



Due to (|211 and ([471) . the matrix \l/ for dynamics relation (jl9|) is now given according to 

/ 1 
* = I 10 

\asA 2 /2 1-8, 
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and u is related to parameter s according to (TT6l) : cr = - — | . Then equivalent mass 

conservation at the order 5 looks like equation (l38l) . We have precisely: 

(48) 

dp dq A 2 At 2 . d 3 q A 4 At 3 <9 4 p 



dt <9x 



1 



12 



a 1 



&r 4 



12 v ~ &r 3 

A 4 At 4 . . , s 2 , d 5 q 

+ — (1 -a) (l + a + 10(1 - 2 a) a 2 ) ' 



120 



dx b 



0(At 5 



Conservation of momentum takes the form: 

— — h « A A At (1 — a) cr 

at :1 



(49) 



dx 



+ C4 



9x 2 
A 4 At 3 8 4 q 

12 ax 4 



+ 



a 4 At 2 a 3 p 



+ 



6 dx 3 
A 6 At 4 <9 5 p 

120 



0(At 5 



with parameters (3 to (5 given by 

C4 
C5 



a(l-a)(l - 6cr^ 
- (1 - a) a (l - 4a - 12 (1 - 2 a) a 2 ) 
a (1 - a) (1 - 4 a - 10 (5 - 9 a) a 2 + 120 (2 - 3 a) a 4 ) . 

When a = -^=, the coefficient (3 of relation (1491) is null. In this case, the lattice Boltzmann 
scheme is formally third-order accurate for the momentum equation. But, as remarked in 
[T6] . the mass conservation (T48l) remains formally second-order accurate, except for the 
(without any practical interest as it leads to a null viscosity) case a = 1. 

• D2Q9 for linearized athermal Navier-Stokes at order four 

The D2Q9 lattice Boltzmann scheme can be used also for simulation of fluid dynamics. 
For the particular case of conservation of mass and momentum, we just replace matrix 
^ of (j4"3i) by the following one, assuming the aim is to simulate an athermal fluid with 
speed of sound \/l/3: 
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We have conservation of mass at fourth-order of accuracy: 



8% 

dy 



A 2 At 2 A 



dq x dq v \ A 4 At 3 . . . 9 



0(At' 
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and conservation of two components of momentum: 
(51) 

dq x A 2 dp A 2 Ai 
dt 3 dx 3 



dx \ dx 



A 4 At 2 



< — 



27 



3(<7 2 + <7?)-l 



d . A 4 At 3 

— An 

dx H 108 



+C: 



d 4 q, 1 



dy J 
d*q x 



+ 07 Ag 2 



10 



22 



+ C 



<9x 2 <9y 2 ' 1513 dx dy 3 ' " 5U * <% 4 
d_(d<h dq y \ 



dx 4 
d 4 q y 



+ C: 



31 



01 



(52) 
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A 4 At 2 



27 



+ cr 7 Ag y 



dg y A 2 <9p . . ^ 
dt ^ dy 3 L dy V cte 

3 ^ 3 + ^ " 1 ) dy- A ? - -108" fl£T + ^ dx^dy 



+V22 



dx 2 dy' 



d^qx . d\ 



dxdy 3 



dy 4 



where the coefficients ( are given by 

C40 = 



(53) 



Vo4 = -cr 3 - a 7 - 12cr 2 (T 7 



12 03 2 + 18 2 5 



+ 6 050 2 — 12 cr 3 04 cr 5 
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0(At 4 



0(At 4 



12 of 



24 3 <7g 7 + 12 4 (7 5 CT 7 
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— 12 03 of — 12 030405 + 12 03 05 o 7 + 12 04 05 07 

r/22 = — 13 03 + 6 a 4 — 10 cr 7 + 18 of 05 — 12of0 7 — 12 3 2 
+ 3Oo 5 0f — 12 03 04 05 + 120 03 05 7 — 60 04 05 7 — 12 3 

77 31 = — 10 03 + 6 04 — 7 7 + 18 u\ 05 — 12 0§ 7 — 12 03 7 
+ I8050 2 + 12 03 04 05 + 84 03 05 7 — 60 04 05 7 + 12 o 3 

r/40 = — 3 7 + 24 05 2 — 12 3 . 



• D3Q19 for linearized Navier-Stokes 

The D3Q19 Lattice Boltzmann scheme is described with details e.g. in J. Tolke et al 
The construction of matrix M that parameterizes the transformation (f41) is presented in 
full detail in relations (1521) to (1571) in Appendix B. The associated matrix \I/ is also of 
order 19 and therefore quite difficult to write on a A4 paper sheet. Due to constitutive 
relations (fl~9l and (1201) . it is easily obtained from the expression of equilibrium moments. 
We have taken for this D3Q19 scheme 
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In order to obtain physical equations at first-order of accuracy with a sound velocity Cq 
given by c = a A the relation 9 = 57 a 2 — 30 must be imposed to obtain correct fluid 
second-order partial differential equations and the parameter (3 remains free. 

• When the number of velocities of the Boltzmann scheme is reduced (up to D2Q9 
scheme typically), it is possible to expand the dispersion equation formally and to derive 
equivalent partial differential equations up to an arbitrary order. We have done the com- 
parison for one-dimensional and bi-dimensional schemes. The process has been extended 
to models with more velocities and various conserved quantities; however the equations 
become very complicated and thus will not be reproduced here. Let us just mention that 
the expressions found are quite similar to those obtained for the previous test cases. 



4 The fourth-order accurate lattice Boltzmann scheme 

• In this section, we precise how to choose particular "quartic" values of relaxation pa- 
rameters in order to increase the accuracy of the scheme. We verify with the help of precise 
numerical experiments for analytical test cases that the numerical precision follows our 
prediction. We focus first on classical thermics at two and three space dimensions. Then 
we propose two numerical experiments for athermal linearized Navier-Stokes equations at 
two and three space dimensions for a nontrivial geometry. 

• The D2Q5 lattice Boltzmann scheme for a thermal problem 

We obtain the order 4 by setting k 40 = and k 2 2 = in relations (TTO and 
respectively. We obtain : 

, rr . a + 4 1 2 + 3a 1 
(55) a 3 = oi- — , 4 



1 — a 12 0i 1 — a 6 0i 

The BGK condition 01 = 03 = 04 leads to 01 = -7= and a = — 4 and thus to a thermal 



/12 

diffusivity equal to 0. Note that the intermediate TRT presented in Ginzburg et al [221123] 
supposes simply 03 = 04. If we insert this constraint inside relations (1551 . we get 

" = 73' = 73 

to enforce fourth-order accuracy. Then the d'Humieres version of lattice Boltzmann 
scheme is mandatory for this improvement of the method with a wide family of admissible 
parameters. In order to study the fourth-order accuracy of the D2Q5 lattice Boltzmann 
scheme for thermal problem, we use three different approaches. The first two consider the 
interior scheme and the third one incorporates boundary conditions. 



• First of all, we study homogeneous plane waves with a "one point computation". In 
that case, we can derive numerically a dispersion equation for scheme (TT5l associated 
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with dH), (fT9l) . (1791) and (1391) . as proposed in [30]. We introduce a wave in the Boltzmann 
scheme, id est f(x, t) = f(k x , k y ) exp(ik x x + ik y y). Then we have f(x, t + At) = 
G f(x, t) with the so-called amplification matrix G (see e.g. Richtmyer and Morton 
P0~] ) obtained without difficulty from matrices M, \l> and S defined respectively in (1791) 
([39|) and 

£ = diag (l, e l k * Ax , e l k v Ax , e~ 2 Ax , e _i fc » Ax ) 

for the D2Q5 scheme displayed in Figure ED (left). Then G = BM^^M. Then if 
Jjj is formally given by relation (flDl and operators ^ and ^ replaced by i k x and 
ifcy respectively, the number z = exp(At^) is an eigenvalue of matrix G at fourth- 
order of accuracy. The numerical experiment (see Figure Q]) confirms the theoretical 
development of the dispersion equation. Note that for situations relaxing to uniform 
state, the eigenvalues that we determine below are negative; however we shall express 
results in terms of positive relaxation rates with adequate sign changes. 

• For inhomogeneous situations, with Nc lattice points (and Nc (J + 1) degrees of 
freedom), one can study the time evolution starting from some initial state. Another 
approach for linear situations considers that the state X(t) that belongs to M. Nc ( J+1 ) 
can be decomposed as a sum of eigenmodes of the operator A defined using the discrete 
evolution scheme: 

(56) X(t + At) = A.X(t). 

The matrix A being of very large size, one can look for some of its eigenmodes using 
for instance the method proposed by Arnoldi [2j. To accelerate the Arnoldi computa- 
tions, following a suggestion by L. Tuckerman [17], we replace the determination of the 
eigenvalues of equation (1561) by the determination of the eigenvalues of 

(57) X(t+ (2£ + 1) At) = A 2e+1 .X(t), 

for some £ G N, using the fact that the lattice Boltzmann scheme is very fast compared to 
the inner "working" of the Arnoldi procedure. Replacing problem (1561) by problem (1571) not 
only increases the splitting between various eigenmodes, but also helps to discriminate 
against the acoustic modes by multiplying the logarithm of the imaginary part of the 
eigenvalues by 2£ + 1. Note that choosing an even number of time steps would bring in 
the "checker-board" type modes. We denote by r num any eigenvalue computed with this 
methodology. 

• We first test this method for "internal" lattice, id est with a periodic N c = N x x N y 
situation and find the same results as those derived from the "one-point" analysis (see 
Figured]) with very good accuracy. For this periodic situation, the eigenmodes are plane 
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Figure 1: Error | ^f^- — 1 1 of D2Q5 scheme for thermic test case, "one point" simulation. 
Different curves correspond to different orientation of the wave- vector with respect to the 
axis, showing the angular dependence of the next order. 



waves for the wave vector k x = ^jf-, k y = 2-71"-^-, where I x and I y are integers. We 
compare the numerical relaxation rates r num (/ x , I y , N x , N y ) to T t h = n(k^. + ky) and show 
in Figure [2] the relative difference between those two quantities (called the "error") for the 
particular values I x = 5 and I y = and N x from 11 to 91. With arbitrarily chosen values 
of the "non-hydrodynamic" s-parameters, we observe second-order convergence. However 
for the quartic s-parameters the convergence is of order four with a large decrease in the 
absolute value of the error. Analogous results are displayed in Figure [2] for D3Q7. 

• We now consider a second case with boundary conditions: exact solution for the 
modes of the Laplace equation in a circle of radius R with homogeneous Dirichlet bound- 
ary conditions. Density is defined with ([ED applied with J = 4 in this particular case. 
Recall that density follows heat equation |j — k Ap = with k = G\ (4 + a) and 
homogeneous boundary conditions at r = R. The solution of this problem is standard 
(see e.g. Landau and Lifchitz [32], Abramowitz and Stegun [TJ or Carslaw and Jaeger 
[6]) and is parameterized by a pair (£, n) of integers. We introduce the n th zero 
of the Bessel function J^. Then a solution with time dependence as exp(— Yt) defines a 
corresponding eigenvalue Y (also denoted as r th in the following) that satisfies 

(58) r = « (|) . 
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Figure 2: Arnoldi test case for periodic thermics, I x — 5, I y — 0. Various parameters for 
lattice Boltzmann schemes D2Q5 and D3Q7. 



• The effect of fourth-order accuracy Boltzmann scheme in computing the eigenfunc- 
tion is spectacular: just compare Figures [3] and QJ Nevertheless, the effect of boundary 
conditions (we use anti-bounce-back with interpolation a la Bouzidi et al [2]) cannot be 
neglected. In Figure El we have compared the error defined by | — 1 1 for two internal 
schemes (with usual and quartic parameters) and two versions (first and second-order) of 
simple numerical boundary conditions introduced by Bouzidi et al [5]. We still observe a 
better numerical precision of the schemes (by two orders of magnitude typically) whereas 
the convergence still remains second-order accurate. We conclude that the effect of bound- 
ary conditions is crucial for the determination of the order of convergence. Nevertheless, 
the choice of quartic parameters gives a higher precision for the lattice Boltzmann scheme. 

• The D3Q7 lattice Boltzmann scheme for a thermal problem 

We obtain the order 4 by setting k 40 o = and k 2 2o = in relations (1451) and (T46l) . We 
obtain : 

1 a + 6 4 + 3a 1 

o ox 1 — a 12 (1 — a) <j\ 

As for D2Q5, the "BGK condition" o\ = 04 = a% leads to a§ = and a = — 6 and thus 
to thermal diffusivity equal to 0. Theoretical modes of the Laplace equation in a sphere 
of radius R with homogeneous Dirichlet boundary conditions are parameterized through 
the n th zero ^™ +1/ / 2 °f semi-integer Bessel function Ji+1/2 and the eigenvalue T is given 
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Figure 3: D2Q5 scheme for thermics inside a circle. Eigenmode n = 4, £ = for 
heat equation with Dirichlet boundary conditions. Second-order accuracy with usual 
parameters for lattice Boltzmann scheme. 
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parameters for lattice Boltzmann scheme. 
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5. Errors for 



by: 
(59) 



Ve+i/2 
R 



e N, n > 1 



• The results of Figures [6] and [7] have been obtained with R = 17.2 and n = 5. The 
theoretical value of the eigenvalue is T = 5 2 n 2 n/R 2 (as for m = 0, the zeros of the 
semi-integer Bessel function are simply nn). We have used parameters Si = 1.26795, 
S4 = 1.2, sq = 1.3 for the usual computations. The quartic parameters have been chosen 
as 

si = 1.26795 , s 4 = s 6 = 0.92820 . 

>From results presented in Figure [8} the conclusion is essentially the same as that ob- 
served for two-dimensional thermics: the results are improved by two orders of magnitude 
typically, but the rate of convergence cannot be rigorously measured or still remains of 
second-order. 

• We also made a parameter study of the location of the boundary condition. We plot 
in Figure [9] the ratio TR 2 /(k,tt 2 ) with T given by relation (1591) . We use Bouzidi et al 
[5] boundary procedure with linear interpolation. The fluctuation due to the boundary 
algorithm is around 0.2 %. The gap between second-order usual computation and new 
fourth-order computation is of the order of 2%. We observe that this gap is one order of 
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Figure 6: D3Q7 lattice Boltzmann scheme for thermics in a sphere. Eigenmode n — 5, 
£ = 1, m = with usual parameters. 




Figure 7: D3Q7 lattice Boltzmann scheme for thermics in a sphere. Eigenmode n — 5, 
£ = 1, m = with quartic parameters. 
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Figure 8: D3Q7 scheme for thermics in a sphere with Dirichlet boundary conditions. 
Eigenmode n = 1, £ = 0. Errors for various parameters for lattice Boltzmann and 
boundary schemes. 
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magnitude larger than the error due to the choice of the boundary condition estimated 
from the fluctuations with the imposed radius. 

• The D2Q9 for linearized athermal Navier-Stokes at fourth-order 

We consider now the linear fluid model obtained by a D2Q9 lattice Boltzmann scheme. 
The equivalent partial differential equations are given at the order 4 by relations ([50 to 
(153|) . The dream would be to enforce high order accuracy. However, this is definitively 
impossible in the framework considered here due to the never null third-order term for 
mass conservation ([501) . Recall notation ((5|) for conservative variables: W = (p, q x , q y )^ 
and write the equivalent equations (l5Tl ) - (l53l) in the synthetic form: 

(60) d t W k + = 0(^ 4 ) • 

h P. <3 

We search for a dissipative mode, id est a mode for linear incompressible Stokes problem 
under the form W(t) = e - rt + l ( k ^ x + k yv) W . Then T is an eigenvalue of the matrix A 
defined by 

A{ = ^2 ^kpq (i k x ) p (i k y ) q . 
h p, <? 

We know (see e.g. Landau and Lifchitz [32]) that for Stokes problem (incompressible 
shear modes), the relation 

(61) r = v (k 2 x + ky) 

is classical. Moreover, as a consequence of (l5B and (|52l 



A 2 

(62) v = —Ata 7 

3 

for a lattice Boltzmann scheme with multiple relaxation times. 

• We propose here to tune the parameters se in such a way that the relation (16D is 
enforced for the modes of (1601 . Precisely, we search si such that 



(63) A m = det 



A - (yAta 7 ) (k* + kl) Id 



0(At 7 ) 



With an elementary formal computation, the third-order term of A m relative to At 
is equal to 



Ay- 3 x 6 / 
(64) > = -^f °7 (kl + ky) ((-1-4 ^ - 8 a 5 a 7 ) + h& 



+ 2(1 - Aa 2 7 - Aa 5 a 7 )k 2 x k 2 y 
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It is then clear that the expression (1641) is identically null for parameters 05 and u 7 
chosen according to 



(65) 



v/3 



05 = -5-, 07 ,, 
3 6 



With this particular choice of parameters, so-called quartic in what follows, the viscosity 
v in relation (1621) has the following particular value: 



(66) 



A 2 At 



v 



0.096225 A At 



Then it is very simple to verify that the determinant A m is null up to terms of order 
seven and relation (1631) is satisfied. 

• As in the particular case of D2Q5 scheme, we have verified with periodic boundary 
conditions that the relaxation rate of a transverse wave is determined with error of order 
six and relative fourth-order precision, as shown in Figure [TU1 The detailed numerical 
convergence plot is very similar to Figure El 

• We have also validated our results for eigenmodes of the Stokes problem inside a 
circle. With the notations introduced previously, the eigenvalues T are given [32] by 



(67) 
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Figure 11: D2Q9 scheme for linear Navier-Stokes. Eigenmode n = 5 £ = 1 for the Stokes 
problem in a circle. 

The result for R = 30.07, £ = 1 and n = 5 is presented in Figure [TT] for the velocity field 
with a mesh included in a square of size 61 x 61. The alternance of directions for the vector 
field is clearly visible on Figure E] and we use around seven meshpoints between two zeros 
of the Bessel function. We have compared with the same mesh the results obtained with 
the lattice Boltzmann scheme with the usual parameters that does not satisfy relation fl65l) 
but such that v = which is very close to f|66l) and quartic parameters. The radial 
profile of the tangential velocity is shown in Figures [12] to [HI The difference is visually 
spectacular. As for the thermics case, we observe that simple boundary conditions (here 
we use those of Bouzidi et al. [5]) prevent fourth-order convergence for the Stokes problem. 
Use of more sophisticated boundary conditions (see Ginzburg and d'Humieres [21]) may 
help to improve the convergence; however for models with limited number of velocities, it 
is not clear whether the choice of s-parameters will be the same for "fourth-order volume" 
and "accurate Poiseuille type boundary conditions". 

• The D3Q19 for linearized athermal Navier-Stokes at fourth-order 

The D3Q19 model is analyzed as was done above for the D2Q9 model. We detail in 
Appendix C the way to enforce the precision of eigenmodes for the Stokes problem. We 
obtain a set of eight equations for the coefficients ex's. These equations have only one 
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Figure 12: D2Q9 scheme for linear Navier-Stokes in a circle. Eigenmode n = 5, t 
for the Stokes problem. Usual parameters. 
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nontrivial family of solutions given according to 



(68) ^ 



energy 


<T 4 = 


1 1 

s 4 2 


«4 


= ad libitum 


stress tensor 


05 = 


1/Vl2 


s 5 


= 3-V3 


energy flux 


010 = 


1/V3 


^10 


= 4^-6 


square of energy 


013 = 


i i 

S13 2 


Sl3 


= ad libitum 


other moments of kinetic energy 


Cr 14 = 




S14 


= 3-V3 


third-order antisymmetric 


Cl6 = 


1/V3 


5l6 


= 4^3-6. 



Note these results are incompatible with BGK hypothesis (all a equal) but are compatible 
with the "two relaxation times" hypothesis which enforces equality of even moments 04 = 
05 = °"i3 = 014, an d of odd moments: a w = a ie . We remark that the relaxation rate 
for energy (linked to the bulk viscosity) is not constrained. Note that the shear viscosity 
v takes the value 1/V108 as in (1661) . As for D2Q9 there is no decoupling at order 3 
of shear and acoustic modes, and thus, at least at the present stage, we make no claim 
concerning possible improvements for the acoustic modes. We will study this question in 
a forthcoming contribution. 

• We have performed the same kind of numerical analysis as for the two-dimensional 
D2Q9 case and find quite similar results. We illustrate our results first with a "one point 
experiment". We introduce numerical wave vectors k close to zero and compute the 
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Figure 15: D3Q19 for "one point" experiment and various directions of the wave vector. 

eigenmodes numerically. The shear mode is close to 4f °s | k | 2 and we plot in Figure fT5l 
the experimental error. With ordinary coefficients, the error is of order 4, whereas with 
the so-called "quartic coefficients", the error is of order 6 and the relative error of order 4. 

• We also illustrate our results for the problem of Stokes modes in a sphere which has 
an analytical solution in terms of Bessel functions. The Stokes problem searches for a 
velocity field u(r, t) with u = on the surface of a sphere of radius R. An analysis similar 
to that for the Stokes problem in a circle, leads to an eigenvalue problem, with solutions 

' /~n \ 2 
V 



3^+1/2 



> 1 



analogous to fl59l) . with C^+1/2 e Q ua l to the n zero of the "semi-integer" Bessel function 
Je+1/2 as defined in Abramowitz and Stegun [T]. Using the Arnoldi technique, we can 
determine a few eigenvalues and verify that they are close to the theoretical formula. We 
present Figures [16] and [T7| an example of a typical result obtained with this framework. 
We find that these eigenvalues have the expected degeneracy 2£+ 1. Note however that, 
the computations being made for a rather small radius R, there are small splittings of the 
degenerate eigenvalues due to the fact that lattice Boltzmann computations have cubic 
symmetry. 



• For a more detailed analysis, we take advantage of the symmetry of the Stokes problem 
and therefore perform computations on an eighth of the sphere, taking proper account 



30 



Frangois Dubois and Pierre Lallemand 




Figure 16: D3Q19 for linear Navier-Stokes in a sphere. Eigenmode n = 3, £ = 1 for 
Stokes problem with Dirichlet boundary conditions. Tangential velocity vector field for a 
plane through the center of the sphere. 




Figure 17: D3Q19 for linear Navier-Stokes in a sphere. Eigenmode n = 3, £ = 1 for Stokes 
problem with Dirichlet boundary conditions. Tangent vector field for a plane orthogonal 
to vector (1, 1, 1). 
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Figure 18: D3Q19 for linear Navier-Stokes in a sphere. First eigenmodes for stationary 
Stokes problem with Dirichlet boundary conditions. 
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of the symmetry with respect to the planes perpendicular to the coordinates x, y, z , 
through the center of the cube (symmetry or anti-symmetry). Using four different com- 
binations of symmetries on the planes, we can determine all the eigenvalues, the other 
combinations leading to the same eigenvalues with only a permutation in the coordinates 
for the eigenmodes. Note that due to the rather high complexity of the Arnoldi prodecure, 
this allows a reduction in computer time of two orders of magnitude. 

• We present in Figure [18] the effect of boundary conditions for a number of values of 
the radius from 29 to 30. We give in Figure fT9l some details for R between 19 and 20 for 
the m = 1, n = 6 mode. There are two sets of data, one for usual s-parameters 

S4 = 1.3, S5 = 1.25 , sio = 1.2, Si3 = 1.4, S14 = 1.25 , sie — 1.3 

and one for the quartic s-parameters given precedingly in (l68l) with 

s 4 = 1.3, si3 = 1.4. 

Similar work has been done for a cube. The results are published in Leriche et ol in [34] . 

5 Conclusion 

• The expansion of equivalent equations that are satisfied by the mean quantities de- 
termined by the lattice Boltzmann method has been described in this contribution and 
explicit formulae given for a few models up to order four in space derivatives. Extending 
either to more complicated models or to higher order derivatives is very simple and does 
not imply new conceptual developments, in particular careful treatment of non commut- 
ing terms that appear in the Chapman-Enskog procedure. The developments imply only 
simple algebraic manipulations that can be performed by a "formal language" program, 
as used here. Note that these developments have a rather high complexity as seen by the 
fact that each order takes roughly ten times as much computer time as the preceding one. 

• With the Taylor expansion method, we can obtain explicit formulae which, then, 
unable us to tune some parameters of the lattice Boltzmann scheme initially proposed 
by d'Humieres in order to capture, up to fourth-order accuracy, shear waves. Of course, 
this extra-precision obtained with a classical scheme is possible only if the viscosity is 
essentially fixed and the expansion done around zero velocity. Even though very few 
situations were studied here, it can be said that tuning the accuracy of the "internal 
code" independently from the method to take care of boundary conditions allows us 
to get useful information concerning these two sources of errors in lattice Boltzmann 
simulations. Future extension of this work will be to try and discriminate between some 
of the numerous proposed ways to deal with boundaries to be able to estimate their 
contributions to errors in comparison to those due to the "internal code". 
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Appendix A. Taylor expansion method 



• We start from relation (TT51) for iteration of the lattice Boltzmann scheme and take 
the momentum of order k. Then 

J 

m k (x, t + At) = Mu fe( x ~ v e At ' *) 

1=0 

j j 

= J2J2 Mm M tp m K x - ve At ' £ ) due to ^ 

1=0 p=0 

J J J 

1=0 p=0 r=0 

due to (HHJ). We have 

j j j 

(69) m k (x, t + At) = Y Yl Yl Mke M ^v ^P r m r (x-v e At,t), < k < J . 

1=0 p=0 r=0 

We expand now momentum m r (x — Vj>At, t) with a Taylor formula of infinite length: 

~°° (At)" - '' 



(70) m r (x-v £ At,t) = En^(~E Maf9a ) ?mr ( 1 ' 

q=0 " ' 

Then due to (IS9l. (jTO]) and (j22l. we have 



q=0 H <x=l 



J J J ^ | 7 | 

(71) m fc (x, t + At) = J2 J2 J2 Mki M tv ^Vr ~j Tj~ ^7 ^"V , < fc < J . 

7 £=0 p=0 r=0 ' ^ ' " 



We can also expand the left hand side of (ffB and we have finally 

00 \+q J J J . | 7 | 

(72) ^_a«m fc = E E E E M « TT ^ ' 0<k<J. 

q=0 7 1=0 p=0 r=0 ' ^ ' ' 



• We consider relation (1721 at order zero relative to time step for a conserved component 
of momentum (id est < k = i < N—l). The left hand side of (1721 is equal to mj + 0(At) 
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and we have 

j j j 

Wi + O(At) = J2J2J2 M ^ M ^v V m r + O(At) 

£=0 p=0 r=0 
J J 

p=0 r=0 
J 

= ^2 ^ir m r + O(At) with < Z < N 

r=0 
./ 



^ 5 ir m r + O(At) due to ([20 



r=0 

= mi + O(At) 

and no information is contained at this first step. Consider now the same development 
for k > N. We pass over some repeated summations: 

m k + 0(At) = M M M' 1 m pr m r + O(At) 

N-l 

= Mke M Zp m i + Yl Mke M tv ^vr ™r + O(At) 

j=0 r>N 
N-l 

= 5k P *W W i + Yl Mk£ M ^V 5 PT i 1 ~ S r) m r + 0( At) 

j=0 r>N 

due to (HI) and (J2DJ 

N-l 

= **i Wj + M ke M h l (l-s p )m p + O(At) 

3=0 

N-l 

= hp (1 - s p ) m p + **i Wj + O(At) 

j=0 

N-l 

= (1 - s fc ) m k + W 3 + °( At ) • 

3=0 



We deduce from the previous calculation the relation (1231) with the expression fl24l) of the 
coefficients . Then we can go now one step further. 

• At first order, relation (f72l becomes 



(73) m k + At-^ + 0(At 2 ) = m* k - At ^ M fc ^ M^" 1 \l/p r M at d a m r + 0(At 2 



dt 

a=l 
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For conserved variables j5]) (id est < k = i < N — 1), we have after dividing by At, 
+ O(At) = - M n M ip *pr M ae d a m r + O(At) 

a=l 
d 

= - ^ Ki V d a m r + O(At) due to (EEBD 

= E A - ( E *PJ ^ + E ^ ^) + ( At ) 
d 

= E A - E + E^t*^)) + °( At ) 

a=l j<N 1>N l 

N-l d J 1 

E E E A - (>• + E **) w + 0(At) • 



j=0 a=l p=0 i>N 

For an index 7 between 1 and d, we define A]j according to the relation (1261) and the 
previous calculation can be written as a conservation law at first-order 

(74) f + = °( A *)> 0<i<JV-l. 

l7l=l 

• We start again from relation (T73l) with nonconservative indices k {k > N): 

dm JV ~ 1 d 

m k = -At + (l-s fc ) m fc + ^kj Wj - At ^ M u M£ V pr M al d a m r + 0(At 2 ) . 

j=0 a=l 





). 










-At 


i* 












At 








+ — 











dWi 

m k = — 1 y/ fc7 - w,- - At — v ki — — 



AtAl k ^ pr d a (^ rj W^ + 0(At 2 ) 

A^V*nW) + 0(At 2 ). 



l7l=l 



We introduce B k - for |7|= 1 according to (1281) and due to previous calculation, relation 
f|23l) can be extended as 

(75) m k = Athl B Ij d i W i + °( At2 ) ' 

o<M<i 

• We generalize the relations (1741) and (TTBD at the order a through a recurrence hy- 
pothesis (l29j) (1301) . In order to treat the left-hand side of relation (1721) . we observe that 
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we have 



d 2 t Wi = - J2 A * l7hl ^7i ^ [p t w^ + 0(At a ) 

1<| 7 |<(7 



1<|5|<<T 1<|e|<<t 

and if we introduce C]f according to (I3T1) and 

C-f = - A ^ A % > 2<| 7 |<a + l, 

I<5|>1, |£|>l,<5+e=7 

we have for the second time derivative a relation quite analogous to fl29lk 

d 2 t Wi + At ' 7 '" 2 C if d i W i = O(Ar) , < % < N - 1 . 

2<| 7 |<a+l 

This relation can be generalized at an arbitrary order according to 

(76) d q t Wi + Athl ~ q C ?f d i W i = 0{At a ) , < i < N - 1 . 

If relation (1761) is true at order g, we have by differentiation with respect to time, 

d q t +1 Wi = - Athl ~ 9 C W d i (w) + o{At a ) 

<?<l7l<°"+9-l 



J2 AtW- q C«/d s ( ^~ l A%d e W>) + O(Af) 

5|<a-+g-l l<|e|<<T 

At' 7 ' -9 " 1 C« +1 ' 7 d 1 (dtWj) + 0(At CT ) 



<3+l<l7l<°-+<? 

and relation (T76l) is satisfied at the order q + 1 with C 9+1 ' 7 given by the recurrence 
relation (1321) . In an analogous way, we have 

(77) d q t m k = Athl ~ 9 D kj d i W i + 0(At <7+1 ) , k > N, 

q<h\<v+q 

with -D^ 7 defined according to (1341) . If the relation (1771) is satisfied at order g, we have 
by differentiation relative to time, 



d q t +1 m k = Yl ^\~ q Dl] d.idtW^) + 0(At CT+1 ) 

= - £ At'" 5 '- 9 ^ A^l- 1 ^.^) + 0(At CT+1 ) 



<?<l<5|<cr+<? 1<|e]<<t 

^ At l7h(?+1) £,9+1.7 fl^. + (At CT+1 ) 



■7+l<lTl< cr+q+1 
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with coefficients D^ 1 ' 1 determined according to the relation (1351) . We observe that for 
the particular value (7 |= a + 1 the coefficient D|t ' 7 is well defined for < q < a . In 
other words, the coefficient is well defined for 1 < q < 1 7 1 . 

• We verify now by induction that the recurrence relations (1291) and (1301) are satisfied. 
It is the case at order 1 as we have shown in (1741) and (1751) . We first consider a label i 
such that < % < N — 1. Then according to (1721) . we have at the order a + 2 : 

Wi + At^ + V ^SfWi + 0(At CT + 2 ) = 
Wi + M u M- 1 V pr TJJi P « 9 & ( E At ' £ ' B rj W J + 0(At CT+2 ) . 

1<|5|<ct+1 ' 0<|£|<CT ' 

We use relation ([761) for the left hand side of previous relation. We get after dividing by 
At 

9=2 y ' <7<|-7|<ct+<?-1 

a .|5| + | e |-l 

E M u ^ P« S^- d 5+£ V^- + 0(At CT+1 ) . 

l<|5|<<r+l,0<|e|<CT ' ' ' 

and the relation (1291) is extended one step further with a coefficient AJ- defined for 
I7I = q + 1 by the recurrence relation (l33l) . For the nonconserved moments (k > N), the 
relation (1721) can be written at the order a + 2 as 



CT+1 At? 



A+|<5| / \ 

+ M M M^ l m pr —-P l5 d 5 [ E &t {£l B e rj d E Wj) +0(A^ 

1<|<5|<(T+1 0<|e|<cr ' 



1 - 

|(5|<<t+1 1 " '" V 0<|e|; 

We use the relation (1771) and we deduce: 

9 =1 9<l7l<o-+<7 



^ Atl 5 l+I £ l 
+ E ""mr - M u M- 1 ^ pr P e5 B £ rj d s+£ Wj + 0(At CT+2 ■ 

1<|<5|<<T+I,0<|e|<cr 



\5\\ 

We set, with 1 7 1 = a + 1 , £: > N , < j < N - 1 



^ V 1 ^-^^-^-^1 " i^UI^- 1 1 n^UI^- „ A i , I I / 



l<q<a+l * l<|<5|<o-+l,0<|e|<o-, <5+e=7 

and the relation (l36l) is established by induction. □ 
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Appendix B. 

Notations for classical lattice Boltzmann schemes 

In order to define precisely our results, the numbering of degrees of freedom must be 
defined and we make this point precise in this Appendix with the help of usual graphics. 
The choice of moments, id est the M matrix (relation (j3J)) is also made explicit. 

i o 2 
• • • 

x- Ax x x+ Ax 

Figure 20: Stencil for the D1Q3 lattice Boltzmann scheme 



• D1Q3 for advective thermics 

Recall first that the D1Q3 lattice Boltzmann scheme (J = 2 in relation (0])) uses three 
neighbours for a given node x: the vertex x itself and the first neighbours located at ±Ax 
from x (see Figure [20j) . We introduce A as in (fTTjl and adopt a labelling for matrix M of 
relation as in Figure [201 

/iii 

(78) M = -A A 

\A 2 /2 A 2 /2 




• D2Q5 for classical thermics 

We have now four (J = 4) nontrivial possible directions for propagation of particles 
(Figure EU left). We adopt for the M matrix of relation (jlj) the following choice: 



(79) 



M 



( 1 





-4 



1 
A 

1 
1 



1 

A 
1 



1 

-A 

1 
1 



1 \ 



-A 
1 

-v 



• D2Q9 for classical thermics 

The lattice Boltzmann model D2Q9 is obtained from the D2Q5 model by adding four 
velocities along the diagonals (Figure EH right). The evaluation of matrix M is entirely 
nontrivial. We refer the reader to |30j, and the reader can also consult our introduction 



Towards higher order lattice Boltzmann schemes 



39 







1 *2 






• 




» 


• 






i 












►o 




Ax 


• 




» 


• 








Ax 







7 4 



Ax 



Ax 



Figure 21: Stencils for D2Q5 and D2Q9 lattice Boltzmann schemes 





Figure 22: Stencils for D3Q7 and D3Q19 lattice Boltzmann schemes 
fT4l. We have: 
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• D3Q7 for pure thermics 

For three-dimensional thermics, one only needs a seven point scheme and can use the 
so-called D3Q7 lattice Boltzmann scheme whose stencil is described in the left part of 
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Figure 1221 The matrix is not very difficult to construct. We follow [3 1 J : 



(81) 



M 



( 1 


1 


1 
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-A 
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1 / 



• D3Q19 for linearized Navier-Stokes 

The D3Q19 Lattice Boltzmann scheme is described with details e.g. in J. Tolke et ol [46] 
and the stencil is presented in Figure [22] (right). The matrix M that parameterizes the 
transformation (jlj) looks like this: 
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Due to the important number of moments, we detail in this sub-section the way the 
previous matrix is obtained. First, velocities t>" for < j < J = 18 and 1 < a < 3 are 
naturally associated with Figure [22l The four first moments p and q a are determined 
according to ([ED and ((71) and the associated elements for matrix M are given in (JEJ) 
and (23). The construction of other moments uses the tensorial nature of the variety of 
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moments that can be constructed, as analyzed by Rubinstein and Luo [41]: scalar fields 
are naturally coupled with one another and similarly for vector fields, and so on. So 
components of kinetic energy are introduced: 

(82) M 4j = 19J2\ | 2 , < j < J . 

a 

The entire set of second-order tensors is completed according to 



M 5j 


= 2(„}) 2 


- K) 2 




2 




M ej 


= - 










M 7j 


= v- ?; 2 

U 3 V 3 ' 


M 8j = 




M 9j = u| v) , 


< j < J 



The three components of heat flux are defined by 



(84) 



v a I 2 



M 10j = hv)Y,\ < r , = E i < i 2 ' = 5u J J2 

a a a 

o<j<J. 

We finally obtain the moments of higher degree: the square of the kinetic energy 

(85) Mi 3j = y (EKI 2 ) 2 ' 0<J<^ 

a 

second-order moments "weighted" by kinetic energy: 

(86) | M 15j = 3 ((v^ - (uj) 2 ) £ I v? | 2 , ° < j < J, 
and third-order anti-symmetric moments: 





[ Miej = 


= *} ((«?)' ~ 


- m 




(87) | 




= -I((^) 2 - 


- M) a ) 






I Ml 8j = 


= t$ ((„})»- 


- (-D 2 ) , 


< J < J 



Then matrix M is orthogonalized from relations P, {82]), {83]), ([84]), (J85J) , (I86l) and 
fl87l) with a Gram-Schmidt classical algorithm: 

Mij = Mij - 9u M tj , i > 4 • 

Ki 

The coefficients ga are computed recursively in order to force orthogonality: 

j 

^2 M ij M kj = for i± k . 

3=0 
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Appendix C. Quartic parameters in three dimensions 

• We use the equivalent equations of lattice Boltzmann scheme D3Q19 obtained pre- 
viously in the following way. We consider the vector of conserved variables ((5j): W = 
(p> Qxj Qy, Q z ) ■ We write the equivalent partial differential equations under the combined 
form: 

(88) d t w k + ^flETO = °( A * 4 )- 

We search dissipative mode solution of ((HSJ) under the form W(t) = e ~ rt+i ( fc * a+ **s' +fc ** ! ) 
Then Y is an eigenvalue of the matrix A defined by 

4 = E A U {ik x y{ih y y{ik z y. 

Pi <?; r 

• We wish to solve this dispersion equation with a high order of accuracy, id est in our 
present case: 

(89) A = det [A - rid] = 0(At 7 ) . 

We impose also that this eigenvalue is double as classical for shear waves in three dimen- 
sions [32]: ^ (det [A — Tld]) ~ 0. The first nontrivial term in powers of At for this 
derivative of the determinant is the term of order 3. Then we force 

(90) -j^ (det [A - rid]) = 0(At 4 ) . 

For Stokes problem (incompressible shear modes) and D3Q19 lattice Boltzmann d'Humieres 
scheme, we have [38] : 

(91) r = v\k\ 2 = ¥-Ata 5 (kl + k 2 y + k 2 z ). 

• We solve the set (1591 (1901) (1911 of equations for all values of the time step At . We 
obtain in this way a set of eight algebraic equations: 

2 CT 5 Oio - 4 erf + 6 cr 5 (X w = 1 

80 of - 32 of <tio + 24 of (7 10 o X6 + 12 o M o X6 of - 8 of - 4 of of 
+12 of of 6 - 12 erf cr 14 cr 10 - 12 a 5 cr 16 cr 14 cr 10 + 6 a 5 a 14 of 

-8 0- 5 Cri 6 + 6(T 5 (T^g CTi4 - (Ti4 (7i 6 + (Ti 4 (Jio + 1 = 

< -48 of oio + 44 a\ af + 2000 a\ o X6 + 95 erf - 16 of o M o i0 

+292 <7i4 (Ti 6 erf + 68 erf cri 4 Oi — 272 of cri 6 o i4 — 1032 of ofg cri 4 
+56 erf o M cr? - 320 of - 1048 of cr 10 er 16 + erf 4 + 60 erf of 6 o^ - 16 o 5 cri 6 o^ 4 
+72 of o? 4 o-io a-ie - 8 o 5 o? 4 <t 10 + 24 of o 14 + 12 of af 4 o? - 248 «r| 
-464 of a-ie a-14 aio + 148 of o xo - 1284 oi 6 of + 4284 o\ % of - 20 o 5 0i 4 = 
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( - 1 + 2 o 5 (Tio - 4 of + 6 o 5 (Ti 6 ) (2 0-5 010 + 2 a u Oi 

— 2 of — 10 cr 5 Oi6 — 2 (T14 (7i6 + 3) = 

96 of a-10 + 24 o| o? - 1920 erf ai 6 + 98 of + 24 a^ a u cr 10 + 350 a M a 16 of 
+34 erf au a w + 264 erf tr 16 o 14 - 1524 of o? 6 o i4 + 12 erf (j 14 af Q 
+240 of - 576 o\ a-10 o-ie + <r? 4 + 102 <*l ^le °"i4 " 20 a 5 cri 6 o-? 4 + 36 o\ a 2 u a w a w 
-4 a 5 aio - 24 of au + 6a 5 2 of 4 cr? + 240 o^ - 216 of o X6 o M oi + 72 of oi 

- 1488 ai6 of + 5688 of 6 o^ - 20 o 5 o M = 
-o 5 + 6 oie of + 2 of 010 - 4 of = 
2 erf O10 - 2oi 6 of + o 5 - o 5 oi 6 oi4 + o 5 Oi 4 oi - 12 of = 
10 o 5 oi 6 0-14 + 2 o 5 o u 010 + Ho 5 -oi4 + 8of-82oi 6 of + 6 of O10 = 0. 

These equations have only one nontrivial family of solutions given by (168|) . 
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